/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

/* wpp.c -- direct evaluation of water polarisation profile. References: RDF algorithms in Allen-Tildesley, 6.2 and Rapaport, 4.3 */

#include "dataStructs.h"

extern const VecR *siteCoords, *siteOrientVec;
extern const Site *site; 
extern const VecR region; 
extern const int nSites, stepCount, limitWpp;
extern const real deltaT, kGB, sigTail, timeNow; 
extern const int halfSizeHistWpp; 
extern real *histWpp, meanWppDeltaZ, meanWppBoxHeight; 
extern int countWpp, *nWatsPerWppSlab;

void ZeroWppHistograms(); // external function declaration 

static void PrintWpp( const int sizeHistWpp )
{
  real distanceFromBilayerCentre;
  int bin;
  FILE *wppFilePtr;
  char wppSnapFile[ 20 ];
  char time_ns[ 9 ]; // string containing timeNow, e.g. "100", in ns

  sprintf( time_ns, "%d", ( int ) ( 0.5 + timeNow *  TIME_ns ) );
  wppSnapFile[ 0 ] = '\0';
  strcat( wppSnapFile, "wpp-" );
  strcat( wppSnapFile, time_ns );
  strcat( wppSnapFile, ".dat" );
  wppFilePtr = fopen( wppSnapFile, "w" );
  
  for ( bin = 0; bin < sizeHistWpp; bin++ ) {
    distanceFromBilayerCentre = ( bin + 0.5 ) * meanWppDeltaZ - meanWppBoxHeight / 2.; 
    fprintf( wppFilePtr, "%7.3f", distanceFromBilayerCentre ); // printing abscissa [nm]
    fprintf( wppFilePtr, "%12.3f\n", histWpp[ bin ] ); // printing orientation (dimensionless)
  }
  fflush( wppFilePtr );
}

void EvalWpp( const int sizeHistWpp, const real wppDeltaZ ) // Compute the rotational order parameter as defined in allen & tildesley section 5.7.3
{
  real slabVolume; // z-distance resolved
  int n, bin; // "bin" is the histogram bin number
  VecR zPlus; 

  VSet( zPlus, 0., 0., 1.); // unit vector along z-axis

  slabVolume = region.x * region.y * wppDeltaZ;   // computing the volume over which the density is worked out
  
  DO_SITE { // computing wpp histogram: sorting sites in bins from 0 to sizeHistWpp - 1 
    if ( WATER_TYPE == site[ n ].type ) {
      //      printf("*** wppDeltaZ = %f\n", wppDeltaZ);
      bin = siteCoords[ n ].z / wppDeltaZ + halfSizeHistWpp; 
      histWpp[ bin ]  += VDot( siteOrientVec[ n ], zPlus ); // cosine between the (unit) water dipole and the (unit) z-axis vector
      ++nWatsPerWppSlab[ bin ];
    }
  }
  
  meanWppDeltaZ += wppDeltaZ;  // wppDeltaZ computed in SingleStep()
  meanWppBoxHeight += region.z; 
    
  ++countWpp;

  if ( limitWpp == countWpp ) {  // if happy with the amount of data collected
    for ( bin = 0; bin < sizeHistWpp; bin++ ) { // normalising to have an averaged value
      if ( 0 != nWatsPerWppSlab[ bin ] ) { // if there *are* water in the current slab (this is to avoid dividing by zero)
	histWpp[ bin ] /= nWatsPerWppSlab[ bin ];
      }
    }
    
    meanWppDeltaZ /= countWpp;
    meanWppBoxHeight /= countWpp; 
    
    PrintWpp( sizeHistWpp );
    countWpp = 0;
    ZeroWppHistograms();
  }
}
